################################################################################################
# REPLICATION CODE FOR "POLITICAL AWARENESS AND THE IDENTITY-TO-POLITICS LINK IN PUBLIC OPINION"
# Anne Marie Green and Paige Tsai
# 12.10.22
################################################################################################

# Load libraries:
library(survey)
library(ggplot2)
library(cowplot)
library(patchwork)
library(tidyverse)
library(ggrepel)
library(dplyr)
library(stargazer)

# Load data files and create weighted versions:
anes <- readRDS("anes.rds")
cces <- readRDS("cces.rds")
w.anes <- svydesign(~1,data=anes,weights=~weight)
w.cces <- svydesign(~1,data=cces,weights=~weight)

################################################################################
# SET SEED, SIMULATION DEFAULTS, GGPLOT THEMES, AND CREATE HELPER FUNCTIONS 
################################################################################

set.seed(1982)

# Simulation defaults for each dataset:
aware.p <- seq(from=0, to=1, by=.01)
svymean(~age+female+married+income+educ+rel.atd+region+race+religion, w.anes, na.rm=T)
svymean(~age+female+married+income+educ+rel.atd+region+race+religion, w.cces, na.rm=T)
a.ctrl <- data.frame(race="White", religion="Mainline Protestant", lgbt=0, union.member=0, veteran=0, female=1, age=47.5, married=1, income=3.17, educ=3.01, rel.atd=2.63, region="South", year="2016", aware.p=aware.p)
c.ctrl <- data.frame(race="White", religion="Mainline Protestant", lgbt=0, union.member=0, veteran=0, female=1, age=47.9, married=1, income=2.81, educ=2.95, rel.atd=2.94, region="South", year="2016", aware.p=aware.p)

# Ggplot themes 
theme_idkn <- function(){
  theme(panel.background=element_rect(fill="gray91"), 
        axis.text=element_text(size=8, color="black") ,
        axis.title=element_text(size=12, color="black"),
        plot.title=element_text(hjust=.5, face="bold"))
}

theme_idkn2 <- function(){
  theme(panel.background=element_rect(fill="gray91"), 
        legend.position="none",
        axis.text=element_text(size=8, color="black") ,
        axis.title=element_text(size=12, color="black"),
        plot.title=element_text(hjust=.5, face="bold"), 
        plot.margin=margin(0.1,2,0.1,0.1, "cm"))
}


# Helper function to calculate and store confidence intervals
ans.ci <- function(x){
  cis <- confint(x)
  ans <- cbind(as.matrix(x), cis)
  ans <- round(ans, digits=4)
  ans <- as.data.frame(ans)
  names(ans) <- c("est", "lo", "hi")
  ans
}
################################################################################
# REGRESSION MODELS FOR MAIN ANALYSES
################################################################################

# CCES data, party identity
reg.c.pid7 <- svyglm(pid7~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+
                       age+married+income+educ+rel.atd+region+year, w.cces)

# CCES data, policy preferences 
reg.c.atts <- svyglm(atts~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+
                       age+married+income+educ+rel.atd+region+year, w.cces)

# ANES data, party identity
reg.a.pid7 <- svyglm(pid7~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+
                       age+married+income+educ+rel.atd+region+year, w.anes)

# ANES data, policy preferences
reg.a.atts <- svyglm(atts~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+
                       age+married+income+educ+rel.atd+region+year, w.anes)

###############################################################################
# AWARENESS AND THE IDENTITY-TO-POLITICS LINK FOR OTHER GROUPS
###############################################################################

# Simulate for different identities across range of awareness
# Create plots for Figure 2

# (a) Policy preferences and religious identities
evs.c.atts.relg <- rbind(ans.ci(predict(reg.c.atts, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "religion", "Catholic"), type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "religion", "Evangelical Protestant"), type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "religion", "Jewish"), type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "religion", "Secular"), type="response")))
evs.c.atts.relg$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.c.atts.relg$aware.p <- aware.p  
gg.c.atts.relg <- ggplot(evs.c.atts.relg, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("(a) Religion") +
  geom_text(data=subset(evs.c.atts.relg, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(0.01,0,0,.015,0), lineheight=.7) + #
  coord_cartesian(clip='off') +
  theme_idkn2()


# (b) Policy preferences and racial/ethnic identities
evs.c.atts.race <- rbind(ans.ci(predict(reg.c.atts, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "race", "Black"), type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "race", "Hispanic"), type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.atts.race$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.atts.race$aware.p <- aware.p  
gg.c.atts.race <- ggplot(evs.c.atts.race, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(b) Race and ethnicity") +
  geom_text(data=subset(evs.c.atts.race, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (c) Policy preferences and gender identities
evs.c.atts.gend <- rbind(ans.ci(predict(reg.c.atts, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "female", 0), type="response")))
evs.c.atts.gend$id.group <- c(rep("Women", length(aware.p)), rep("Men", length(aware.p)))
evs.c.atts.gend$aware.p <- aware.p  
gg.c.atts.gend <- ggplot(evs.c.atts.gend, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(c) Gender") +
  geom_text(data=subset(evs.c.atts.gend, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (d) Policy preferences and union membership
evs.c.atts.unin <- rbind(ans.ci(predict(reg.c.atts, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "union.member", 1), type="response")))
evs.c.atts.unin$id.group <- c(rep("Non-\nmember", length(aware.p)), rep("Union\nmember", length(aware.p)))
evs.c.atts.unin$aware.p <- aware.p  
gg.c.atts.unin <- ggplot(evs.c.atts.unin, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(d) Union membership") +
  geom_text(data=subset(evs.c.atts.unin, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(-.005,.005), lineheight=.7) + # 
  coord_cartesian(clip='off') +
  theme_idkn2()


# (e) Policy preferences and veteran status
evs.c.atts.vets <- rbind(ans.ci(predict(reg.c.atts, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.atts, replace(c.ctrl, "veteran", 1), type="response")))
evs.c.atts.vets$id.group <- c(rep("Non-\nveteran", length(aware.p)), rep("Veteran", length(aware.p)))
evs.c.atts.vets$aware.p <- aware.p  
gg.c.atts.vets <- ggplot(evs.c.atts.vets, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(e) Veteran status") +
  geom_text(data=subset(evs.c.atts.vets, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (f) Party identity and religious identities
evs.c.pid7.relg <- rbind(ans.ci(predict(reg.c.pid7, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "religion", "Catholic"), type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "religion", "Evangelical Protestant"), type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "religion", "Jewish"), type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "religion", "Secular"), type="response")))
evs.c.pid7.relg$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.c.pid7.relg$aware.p <- aware.p  
gg.c.pid7.relg <- ggplot(evs.c.pid7.relg, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted party identity") + ggtitle("(f) Religion") +
  geom_text(data=subset(evs.c.pid7.relg, aware.p==1), aes(label=id.group), hjust=0, nudge_x=0.02, nudge_y=c(0.015,-.005,0,0,0), lineheight=.7) + 
  coord_cartesian(clip='off') +
  theme_idkn2()


# (g) Party identity and racial/ethnic identities
evs.c.pid7.race <- rbind(ans.ci(predict(reg.c.pid7, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "race", "Black"), type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "race", "Hispanic"), type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.pid7.race$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.pid7.race$aware.p <- aware.p  
gg.c.pid7.race <- ggplot(evs.c.pid7.race, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(g) Race and ethnicity") +
  geom_text(data=subset(evs.c.pid7.race, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# (h) Party identity and gender identities
evs.c.pid7.gend <- rbind(ans.ci(predict(reg.c.pid7, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "female", 0), type="response")))
evs.c.pid7.gend$id.group <- c(rep("Women", length(aware.p)), rep("Men", length(aware.p)))
evs.c.pid7.gend$aware.p <- aware.p  
gg.c.pid7.gend <- ggplot(evs.c.pid7.gend, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(h) Gender") +
  geom_text(data=subset(evs.c.pid7.gend, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (i) Party identity and union membership
evs.c.pid7.unin <- rbind(ans.ci(predict(reg.c.pid7, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "union.member", 1), type="response")))
evs.c.pid7.unin$id.group <- c(rep("Non-\nmember", length(aware.p)), rep("Union\nmember", length(aware.p)))
evs.c.pid7.unin$aware.p <- aware.p  
gg.c.pid7.unin <- ggplot(evs.c.pid7.unin, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(i) Union membership") +
  geom_text(data=subset(evs.c.pid7.unin, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (j) Party identity and veteran status
evs.c.pid7.vets <- rbind(ans.ci(predict(reg.c.pid7, c.ctrl, type="response")),
                         ans.ci(predict(reg.c.pid7, replace(c.ctrl, "veteran", 1), type="response")))
evs.c.pid7.vets$id.group <- c(rep("Non-\nveteran", length(aware.p)), rep("Veteran", length(aware.p)))
evs.c.pid7.vets$aware.p <- aware.p  
gg.c.pid7.vets <- ggplot(evs.c.pid7.vets, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("(j) Veteran status") +
  geom_text(data=subset(evs.c.pid7.vets, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# Create Figure 2
f2_row1 <- gg.c.atts.relg + gg.c.atts.race + gg.c.atts.gend + gg.c.atts.unin + gg.c.atts.vets + plot_layout(ncol=5)
f2_row2 <- gg.c.pid7.relg + gg.c.pid7.race + gg.c.pid7.gend + gg.c.pid7.unin + gg.c.pid7.vets + plot_layout(ncol=5)
f2_title1 <- ggdraw() + draw_label("Predicted policy preferences", fontface='bold', size=16, hjust=.5)
f2_title2 <- ggdraw() + draw_label("Predicted party identity", fontface='bold', size=16, hjust=.5)

f2_title1 / f2_row1 / plot_spacer() / f2_title2 / f2_row2 + plot_layout(heights=c(.15, 1,.05,.15,1))

############################################
#
#COMPARING CCES AND ANES DATASETS &
#COMPARING 2016 AND 2018 RESULTS
#
###########################################
#anes_2016
anes_16 <- anes %>% filter(year == 2016)
w.anes_16 <- svydesign(~1,data=anes_16,weights=~weight)

#cces_2016
cces_16 <- cces %>% filter(year == 2016)
w.cces_16 <- svydesign(~1,data=cces_16,weights=~weight)

#cces_2018
cces_18 <- cces %>% filter(year == 2018)
w.cces_18 <- svydesign(~1,data=cces_18,weights=~weight)

################################################################################
# SET SEED, SIMULATION DEFAULTS, GGPLOT THEMES, AND CREATE HELPER FUNCTIONS 
################################################################################

set.seed(1982)

# Simulation defaults for each dataset:
aware.p <- seq(from=0, to=1, by=.01)
svymean(~age+female+married+income+educ+rel.atd+region+race+religion, w.anes_16, na.rm=T)
svymean(~age+female+married+income+educ+rel.atd+region+race+religion, w.cces_16, na.rm=T)
a.ctrl.16 <-
  data.frame(
    race = "White",
    religion = "Mainline Protestant",
    lgbt = 0,
    union.member = 0,
    veteran = 0,
    female = 1,
    age = 47.5,
    married = 1,
    income = 3.17,
    educ = 3.01,
    rel.atd = 2.63,
    region = "South",
    aware.p = aware.p
  )
c.ctrl.16 <-
  data.frame(
    race = "White",
    religion = "Mainline Protestant",
    lgbt = 0,
    union.member = 0,
    veteran = 0,
    female = 1,
    age = 47.9,
    married = 1,
    income = 2.81,
    educ = 2.95,
    rel.atd = 2.94,
    region = "South",
    aware.p = aware.p
  )

# Ggplot themes 
theme_idkn <- function(){
  theme(panel.background=element_rect(fill="gray91"), 
        axis.text=element_text(size=8, color="black") ,
        axis.title=element_text(size=12, color="black"),
        plot.title=element_text(hjust=.5, face="bold"))
}

theme_idkn2 <- function(){
  theme(panel.background=element_rect(fill="gray91"), 
        legend.position="none",
        axis.text=element_text(size=8, color="black") ,
        axis.title=element_text(size=12, color="black"),
        plot.title=element_text(hjust=.5, face="bold"), 
        plot.margin=margin(0.1,2,0.1,0.1, "cm"))
}


# Helper function to calculate and store confidence intervals
ans.ci <- function(x){
  cis <- confint(x)
  ans <- cbind(as.matrix(x), cis)
  ans <- round(ans, digits=4)
  ans <- as.data.frame(ans)
  names(ans) <- c("est", "lo", "hi")
  ans
}


################################################################################
# DATASET COMPARISON
################################################################################

# CCES data, party identity -- 2016 ONLY
reg.c.pid7_16 <- svyglm(pid7~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.cces_16)

# CCES data, policy preferences -- 2016 ONLY
reg.c.atts_16 <- svyglm(atts~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.cces_16)

# ANES data, party identity -- 2016 ONLY
reg.a.pid7_16 <- svyglm(pid7~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.anes_16)

# ANES data, policy preferences -- 2016 ONLY
reg.a.atts_16 <- svyglm(atts~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.anes_16)


# CCES data, party identity -- 2018 ONLY
reg.c.pid7_18 <- svyglm(pid7~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.cces_18)

# CCES data, policy preferences -- 2018 ONLY
reg.c.atts_18 <- svyglm(atts~
                          race*aware.p+
                          religion*aware.p+
                          female*aware.p+
                          lgbt*aware.p+
                          union.member*aware.p+
                          veteran*aware.p+
                          age+married+income+educ+rel.atd+region, w.cces_18)

# CCES data, policy preferences -- both years, interaction
reg.c.atts <- svyglm(atts~
                       race*aware.p*year+
                       religion*aware.p*year+
                       female*aware.p*year+
                       lgbt*aware.p*year+
                       union.member*aware.p*year+
                       veteran*aware.p*year+
                       age+married+income+educ+rel.atd+region+year, w.cces)

# CCES data, party identity -- both years, interaction
reg.c.pid7 <- svyglm(pid7~
                       race*aware.p*year+
                       religion*aware.p*year+
                       female*aware.p*year+
                       lgbt*aware.p*year+
                       union.member*aware.p*year+
                       veteran*aware.p*year+
                       age+married+income+educ+rel.atd+region+year, w.cces)

################################################################################
# TABLES 4-6
################################################################################
#Compare the coefficients across ANES and CCES dataset
#to see if they are significantly different
#policy preferences
pvals = NULL
for (i in 1:length(reg.c.atts_16$coefficients)){
  top = (reg.c.atts_16$coefficients[i] - reg.a.atts_16$coefficients[i])
  bottom = sqrt(((summary(reg.c.atts_16)$coefficients[i, 2])^2)+((summary(reg.a.atts_16)$coefficients[i, 2]))^2)
  z = top/bottom
  p = 2*pnorm(abs(z),lower.tail = FALSE)
  pvals = c(pvals, p)
}  

pvals = round(pvals, digits = 3)  
pvals

## FOR PARTY ID
pvals_pid = NULL
for (i in 1:length(reg.c.pid7_16$coefficients)){
  top = (reg.c.pid7_16$coefficients[i] - reg.a.pid7_16$coefficients[i])
  bottom = sqrt(((summary(reg.c.pid7_16)$coefficients[i, 2])^2)+((summary(reg.a.pid7_16)$coefficients[i, 2]))^2)
  z = top/bottom
  p = 2*pnorm(abs(z),lower.tail = FALSE)
  pvals_pid = c(pvals_pid, p)
}  

pvals_pid = round(pvals_pid, digits = 3)  
pvals_pid
pvals

# Main models predicting policy preferences and party identity -- 2016 ONLY
#Policy Preferences
stargazer(list(reg.c.atts_16,reg.a.atts_16),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          dep.var.labels.include = F,
          column.labels = c("Policy (CCES '16)", "Policy (ANES '16)"),
          type = "latex"
)

#Party ID 2016 ONLY
stargazer(
  list(reg.c.pid7_16, reg.a.pid7_16),
  single.row = T,
  stars = c(.001, .01, .05, .1),
  fontsize = "small",
  dep.var.labels.include = F,
  column.labels = c("Party (CCES '16)", "Party (ANES '16)"),
  type = "latex"
  
)

#Main models predicting policy preferences and party identity -- 2016 vs 2018
#Policy Preferences
stargazer(list(reg.c.atts_16,reg.c.atts_18),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          dep.var.labels.include = F,
          column.labels = c("Policy (CCES '16)", "Policy (CCES '18)"),
          type = "latex")

## FOR 2018 Policy Prefs
pvals_atts = NULL
for (i in 1:length(reg.c.atts_16$coefficients)){
  top = (reg.c.atts_18$coefficients[i] - reg.c.atts_16$coefficients[i])
  bottom = sqrt(((summary(reg.c.atts_18)$coefficients[i, 2])^2)+((summary(reg.c.atts_16)$coefficients[i, 2]))^2)
  z = top/bottom
  p = 2*pnorm(abs(z),lower.tail = FALSE)
  pvals_atts = c(pvals_atts, p)
}  

pvals_atts = round(pvals_atts, digits = 3)  
pvals_atts

#Policy Preferences
stargazer(
  list(reg.c.pid7_16, reg.c.pid7_18),
  single.row = T,
  stars = c(.001, .01, .05, .1),
  fontsize = "small",
  dep.var.labels.include = F,
  column.labels = c("Party ID (CCES '16)", "Party ID (CCES '18)"),
  type = "latex"
)

## FOR 2018 PARTY ID
pvals_pid = NULL
for (i in 1:length(reg.c.pid7_16$coefficients)){
  top = (reg.c.pid7_18$coefficients[i] - reg.c.pid7_16$coefficients[i])
  bottom = sqrt(((summary(reg.c.pid7_18)$coefficients[i, 2])^2)+((summary(reg.c.pid7_16)$coefficients[i, 2]))^2)
  z = top/bottom
  p = 2*pnorm(abs(z),lower.tail = FALSE)
  pvals_pid = c(pvals_pid, p)
}  

pvals_pid = round(pvals_pid, digits = 3)  
pvals_pid



#INTERACTED MODEL
#cces 2018 data with interaction
stargazer(list(reg.c.atts,reg.c.pid7),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          dep.var.labels.include = F,
          column.labels = c("Policy (CCES '16)", "Policy (CCES '18)"),
          type = "latex")

###############################################################################
# AWARENESS AND THE IDENTITY-TO-POLITICS LINK FOR OTHER GROUPS
###############################################################################

# (a) Policy preferences and religious identities - CCES 2016
evs.c.atts.relg_16 <- rbind(ans.ci(predict(reg.c.atts_16, c.ctrl.16, type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "religion", "Catholic"), type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "religion", "Evangelical Protestant"), type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "religion", "Jewish"), type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "religion", "Secular"), type="response")))
evs.c.atts.relg_16$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.c.atts.relg_16$aware.p <- aware.p  
gg.c.atts.relg_16 <- ggplot(evs.c.atts.relg_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.25,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Religion") +
  geom_text(data=subset(evs.c.atts.relg_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(0.01,0,0,.015,0), lineheight=.7) + #
  coord_cartesian(clip='off') +
  theme_idkn2()


# (a) Policy preferences and religious identities - ANES 2016
evs.a.atts.relg_16 <- rbind(ans.ci(predict(reg.a.atts_16, a.ctrl.16, type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "religion", "Catholic"), type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "religion", "Evangelical Protestant"), type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "religion", "Jewish"), type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "religion", "Secular"), type="response")))
evs.a.atts.relg_16$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.a.atts.relg_16$aware.p <- aware.p  
gg.a.atts.relg_16 <- ggplot(evs.a.atts.relg_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.25,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Religion") +
  geom_text(data=subset(evs.a.atts.relg_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(0.01,0,0,.015,0), lineheight=.7) + #
  coord_cartesian(clip='off') +
  theme_idkn2()

# (b) Policy preferences and racial/ethnic identities - CCES 2016
evs.c.atts.race_16 <- rbind(ans.ci(predict(reg.c.atts_16, c.ctrl.16, type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "race", "Black"), type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "race", "Hispanic"), type="response")),
                            ans.ci(predict(reg.c.atts_16, replace(c.ctrl.16, "race", "Asian"), type="response")))
evs.c.atts.race_16$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.atts.race_16$aware.p <- aware.p  
gg.c.atts.race_16 <- ggplot(evs.c.atts.race_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.2,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Policy Preferences") + ggtitle("Race and ethnicity") +
  geom_text(data=subset(evs.c.atts.race_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (b) Policy preferences and racial/ethnic identities - ANES 2016
evs.a.atts.race_16 <- rbind(ans.ci(predict(reg.a.atts_16, a.ctrl.16, type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "race", "Black"), type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "race", "Hispanic"), type="response")),
                            ans.ci(predict(reg.a.atts_16, replace(a.ctrl.16, "race", "Asian"), type="response")))
evs.a.atts.race_16$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.a.atts.race_16$aware.p <- aware.p  
gg.a.atts.race_16 <- ggplot(evs.a.atts.race_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.2,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Policy Preferences") + ggtitle("Race and ethnicity") +
  geom_text(data=subset(evs.a.atts.race_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()



## Party IDs

# (c) Party IDs and religious identities - CCES 2016
evs.c.pid7.relg_16 <- rbind(ans.ci(predict(reg.c.pid7_16, c.ctrl.16, type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "religion", "Catholic"), type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "religion", "Evangelical Protestant"), type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "religion", "Jewish"), type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "religion", "Secular"), type="response")))
evs.c.pid7.relg_16$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.c.pid7.relg_16$aware.p <- aware.p  
gg.c.pid7.relg_16 <- ggplot(evs.c.pid7.relg_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.25,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party ID") + ggtitle("Religion") +
  geom_text(data=subset(evs.c.pid7.relg_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(0.01,0,0,.015,0), lineheight=.7) + #
  coord_cartesian(clip='off') +
  theme_idkn2()


# (a) Policy preferences and religious identities - ANES 2016
evs.a.pid7.relg_16 <- rbind(ans.ci(predict(reg.a.pid7_16, a.ctrl.16, type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "religion", "Catholic"), type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "religion", "Evangelical Protestant"), type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "religion", "Jewish"), type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "religion", "Secular"), type="response")))
evs.a.pid7.relg_16$id.group <- c(rep("Mainline\nProtestant", length(aware.p)), rep("Catholic", length(aware.p)), rep("Evangelical\nProtestant", length(aware.p)), rep("Jewish", length(aware.p)), rep("Secular", length(aware.p)))  
evs.a.pid7.relg_16$aware.p <- aware.p  
gg.a.pid7.relg_16 <- ggplot(evs.a.pid7.relg_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.25,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party ID") + ggtitle("Religion") +
  geom_text(data=subset(evs.a.pid7.relg_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, nudge_y=c(0.01,0,0,.015,0), lineheight=.7) + #
  coord_cartesian(clip='off') +
  theme_idkn2()


# (d) Party ID and racial/ethnic identities - CCES 2016
evs.c.pid7.race_16 <- rbind(ans.ci(predict(reg.c.pid7_16, c.ctrl.16, type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "race", "Black"), type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "race", "Hispanic"), type="response")),
                            ans.ci(predict(reg.c.pid7_16, replace(c.ctrl.16, "race", "Asian"), type="response")))
evs.c.pid7.race_16$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.pid7.race_16$aware.p <- aware.p  
gg.c.pid7.race_16 <- ggplot(evs.c.pid7.race_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.2,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party ID") + ggtitle("Race and ethnicity") +
  geom_text(data=subset(evs.c.pid7.race_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (d) Party ID and racial/ethnic identities - ANES 2016
evs.a.pid7.race_16 <- rbind(ans.ci(predict(reg.a.pid7_16, a.ctrl.16, type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "race", "Black"), type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "race", "Hispanic"), type="response")),
                            ans.ci(predict(reg.a.pid7_16, replace(a.ctrl.16, "race", "Asian"), type="response")))
evs.a.pid7.race_16$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.a.pid7.race_16$aware.p <- aware.p  
gg.a.pid7.race_16 <- ggplot(evs.a.pid7.race_16, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(.2,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party ID") + ggtitle("Race and ethnicity") +
  geom_text(data=subset(evs.a.pid7.race_16, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()

f2_row1 <- gg.c.pid7.relg_16 + gg.c.atts.relg_16 + gg.c.pid7.race_16 + gg.c.atts.race_16 +  plot_layout(ncol = 4)
f2_row2 <- gg.a.pid7.relg_16 + gg.a.atts.relg_16 + gg.a.pid7.race_16 + gg.a.atts.race_16 + plot_layout(ncol = 4)
f2_title1 <- ggdraw() + draw_label("CCES 2016", fontface='bold', size=16, hjust=.5)
f2_title2 <- ggdraw() + draw_label("ANES 2016", fontface='bold', size=16, hjust=.5)
f2_title1 / f2_row1 / plot_spacer() / f2_title2 / f2_row2 + plot_layout(heights=c(.15, 1,.05,.15,1))

################################################################################
# DESCRIPTIVE STATISTICS FOR DATASET COMPARISON
################################################################################
# Additional libraries for appendix material
library(gtsummary)
library(texreg)
# Descriptive statistics for CCES data, for Tables A2 and A4
ta2 <- tbl_svysummary(w.cces_16, 
                      type=list(pid7~'continuous', aware.p~'continuous', lgbt='continuous', female='continuous', union.member='continuous', veteran='continuous', married='continuous', income='continuous', educ='continuous', rel.atd='continuous', ideo='continuous', vote='categorical'),
                      statistic=list(all_continuous() ~ "{min} {max} {mean} {sd}"),
                      digits=list(age~2),
                      missing="no")

ta2
ta3 <- tbl_svysummary(w.cces_16, 
                      type=list(lgbt='categorical', female='categorical', union.member='categorical', veteran='categorical'),
                      statistic=list(all_continuous() ~ "{min} {max} {mean} {sd}",
                                     all_categorical() ~"{n} ({p}%) {N}"),
                      missing="no")

ta3

as_kable_extra(ta3, format = "latex")
# Descriptive statistics for ANES data, for Tables A3 and A4
ta4 <- tbl_svysummary(w.anes_16, 
                      type=list(pid7~'continuous', aware.p~'continuous', lgbt='continuous', female='continuous', union.member='continuous', veteran='continuous', married='continuous', income='continuous', educ='continuous', rel.atd='continuous', ideo='continuous', vote='categorical', lf.race='continuous', id.race='continuous', id.relg='continuous'),
                      statistic=list(all_continuous() ~ "{min} {max} {mean} {sd}"),
                      digits=list(age~2),
                      missing="no")
ta4

ta5 <-tbl_svysummary(w.anes_16, 
                     type=list(lgbt='categorical', female='categorical', union.member='categorical', veteran='categorical'),
                     statistic=list(all_continuous() ~ "{min} {max} {mean} {sd}",
                                    all_categorical() ~"{n} ({p}%) {N}"),
                     missing="no")

ta3 
ta5

##Tables
as_kable_extra(ta2, format = "latex")

tbl_svysummary(w.cces, 
               type=list(lgbt='categorical', female='categorical', union.member='categorical', veteran='categorical'),
               statistic=list(all_continuous() ~ "{min} {max} {mean} {sd}",
                              all_categorical() ~"{n} ({p}%) {N}"),
               missing="no")


#######################################################
#RACE AND GENDER INTERACTION
#######################################################
# CCES data, party identity ** interaction race * gender
reg.c.pid8 <- svyglm(pid7~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       race*female*aware.p + 
                       veteran*aware.p+
                       age+married+income+educ+rel.atd+region+year, w.cces)

# CCES data, policy preferences ** interaction race and gender
reg.c.atts_2 <- svyglm(atts~
                         race*aware.p+
                         religion*aware.p+
                         female*aware.p+
                         lgbt*aware.p+
                         union.member*aware.p+
                         veteran*aware.p+
                         race*female*aware.p+
                         age+married+income+educ+rel.atd+region+year, w.cces)

# ANES data, party identity ** interaction race and gender
reg.a.pid8 <- svyglm(pid7~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+
                       race*female*aware.p + 
                       age+married+income+educ+rel.atd+region+year, w.anes)

# ANES data, policy preferences ** interaction race*gender
reg.a.atts_2 <- svyglm(atts~
                         race*aware.p+
                         religion*aware.p+
                         female*aware.p+
                         lgbt*aware.p+
                         union.member*aware.p+
                         veteran*aware.p+
                         race*female*aware.p+
                         age+married+income+educ+rel.atd+region+year, w.anes)

## Race and Gender ## 

# (female) Party identity and racial/ethnic identities
evs.c.pid8.race_female <- rbind(ans.ci(predict(reg.c.pid8, c.ctrl, type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Black"), type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Hispanic"), type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.pid8.race_female$id.group <- c(rep("White", length(aware.p)), 
                                     rep("Black", length(aware.p)), 
                                     rep("Hispanic", length(aware.p)), 
                                     rep("Asian", length(aware.p)))  
evs.c.pid8.race_female$aware.p <- aware.p  
gg.c.pid7.race_female <- ggplot(evs.c.pid8.race_female, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and ethnicity (Female)") +
  geom_text(data=subset(evs.c.pid8.race_female, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (male) Party identity, Race & Gender
c.ctrl_male = replace(c.ctrl, "female", 0)
evs.c.pid8.race_male <- rbind(ans.ci(predict(reg.c.pid8, c.ctrl_male, type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Black"), type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Hispanic"), type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Asian"), type="response")))
evs.c.pid8.race_male$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.pid8.race_male$aware.p <- aware.p  
gg.c.pid8.race_gender_male <- ggplot(evs.c.pid8.race_male, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("Race and ethnicity (Male)") +
  geom_text(data=subset(evs.c.pid8.race_male, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()


interaction_2 = gg.c.pid7.race_female + gg.c.pid8.race_gender_male + plot_layout(nrow = 1)
interaction_2


# (female) Policy preferences and racial/ethnic identities
evs.c.atts.race_female <- rbind(ans.ci(predict(reg.c.atts_2, c.ctrl, type="response")),
                                ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Black"), type="response")),
                                ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Hispanic"), type="response")),
                                ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.atts.race_female$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.atts.race_female$aware.p <- aware.p  
gg.c.atts.race_female <- ggplot(evs.c.atts.race_female, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Policy Preferences") + ggtitle("Race and Gender (Female)") +
  geom_text(data=subset(evs.c.atts.race_female, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# (male) Policy preferences and racial/ethnic identities
c.ctrl_male_2 = replace(c.ctrl, "female", 0)
evs.c.atts.race_male <- rbind(ans.ci(predict(reg.c.atts_2, c.ctrl, type="response")),
                              ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Black"), type="response")),
                              ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Hispanic"), type="response")),
                              ans.ci(predict(reg.c.atts_2, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.atts.race_male$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.atts.race_male$aware.p <- aware.p  
gg.c.atts.race_male <- ggplot(evs.c.atts.race_male, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("Race and Gender (Male)") +
  geom_text(data=subset(evs.c.atts.race_female, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.7) +
  coord_cartesian(clip='off') +
  theme_idkn2()

## CREATE FIGURE ##
interaction_1 = gg.c.atts.race_female + gg.c.atts.race_male + plot_layout(nrow = 1)

## Race and Gender ## 

# (female) Party identity and racial/ethnic identities
evs.c.pid8.race_female <- rbind(ans.ci(predict(reg.c.pid8, c.ctrl, type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Black"), type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Hispanic"), type="response")),
                                ans.ci(predict(reg.c.pid8, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.pid8.race_female$id.group <- c(rep("White", length(aware.p)), 
                                     rep("Black", length(aware.p)), 
                                     rep("Hispanic", length(aware.p)), 
                                     rep("Asian", length(aware.p)))  
evs.c.pid8.race_female$aware.p <- aware.p  
gg.c.pid7.race_female <- ggplot(evs.c.pid8.race_female, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and ethnicity (Female)") +
  geom_text(data=subset(evs.c.pid8.race_female, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()


# (male) Party identity, Race & Gender
c.ctrl_male = replace(c.ctrl, "female", 0)
evs.c.pid8.race_male <- rbind(ans.ci(predict(reg.c.pid8, c.ctrl_male, type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Black"), type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Hispanic"), type="response")),
                              ans.ci(predict(reg.c.pid8, replace(c.ctrl_male, "race", "Asian"), type="response")))
evs.c.pid8.race_male$id.group <- c(rep("White", length(aware.p)), rep("Black", length(aware.p)), rep("Hispanic", length(aware.p)), rep("Asian", length(aware.p)))  
evs.c.pid8.race_male$aware.p <- aware.p  
gg.c.pid8.race_gender_male <- ggplot(evs.c.pid8.race_male, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("") + ggtitle("Race and ethnicity (Male)") +
  geom_text(data=subset(evs.c.pid8.race_male, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

f2_row1 <- gg.c.atts.race_female + gg.c.atts.race_male + plot_layout(nrow = 1)
f2_row2 <- gg.c.pid7.race_female + gg.c.pid8.race_gender_male + plot_layout(nrow = 1)
f2_title1 <- ggdraw() 
f2_title2 <- ggdraw() 
race_gender = f2_title1 / f2_row1 / plot_spacer() / f2_title2 / f2_row2 + plot_layout(heights=c(.3, 1,.05,.3,1))
race_gender

# Table with main models predicting policy preferences and party identity, race*gender interaction
stargazer(list(reg.c.atts_2, reg.c.pid8),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          dep.var.labels.include = F,
          column.labels = c("Policy (CCES)", "Party (CCES)"),
          type = "latex")
          
#######################################################
# RACE AND REGION INTERACTION 
#######################################################
## Race and Region ##
# CCES data, party identity ** interaction race and region
reg.c.pid7_race.reg <- svyglm(pid7~
                                race*aware.p+
                                religion*aware.p+
                                female*aware.p+
                                lgbt*aware.p+
                                union.member*aware.p+
                                race*region*aware.p + 
                                veteran*aware.p+
                                age+married+income+educ+rel.atd+region+year, w.cces)

# CCES data, policy preferences ** interaction race and region
reg.c.atts_race.reg <- svyglm(atts~
                                race*aware.p+
                                religion*aware.p+
                                female*aware.p+
                                lgbt*aware.p+
                                union.member*aware.p+
                                veteran*aware.p+
                                race*region*aware.p+
                                age+married+income+educ+rel.atd+region+year, w.cces)

# ANES data, party identity ** interaction race and region
reg.a.pid7_race.reg <- svyglm(pid7~
                                race*aware.p+
                                religion*aware.p+
                                female*aware.p+
                                lgbt*aware.p+
                                union.member*aware.p+
                                veteran*aware.p+
                                race*region*aware.p+
                                age+married+income+educ+rel.atd+region+year, w.anes)

# ANES data, policy preferences ** interaction race*region
reg.a.atts_race.reg<- svyglm(atts~
                               race*aware.p+
                               religion*aware.p+
                               female*aware.p+
                               lgbt*aware.p+
                               union.member*aware.p+
                               veteran*aware.p+
                               race*region*aware.p+
                               age+married+income+educ+rel.atd+region+year, w.anes)

######### Race and Region #########

# CCES Party identity and interaction with Race and Region (Northeast)
c.ctrl_north = replace(c.ctrl, "region", "Northeast")
evs.c.pid7.race_north <- rbind(ans.ci(predict(reg.c.pid7_race.reg, c.ctrl_north, type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_north, "race", "Black"), type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_north, "race", "Hispanic"), type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_north, "race", "Asian"), type="response")))
evs.c.pid7.race_north$id.group <- c(rep("White", length(aware.p)), 
                                    rep("Black", length(aware.p)), 
                                    rep("Hispanic", length(aware.p)), 
                                    rep("Asian", length(aware.p)))  
evs.c.pid7.race_north$aware.p <- aware.p  
gg.c.pid7.race_north <- ggplot(evs.c.pid7.race_north, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and Region (Northeast)") +
  geom_text(data=subset(evs.c.pid7.race_north, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (South)
evs.c.pid7.race_south <- rbind(ans.ci(predict(reg.c.pid7_race.reg, c.ctrl, type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl, "race", "Black"), type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl, "race", "Hispanic"), type="response")),
                               ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.pid7.race_south$id.group <- c(rep("White", length(aware.p)), 
                                    rep("Black", length(aware.p)), 
                                    rep("Hispanic", length(aware.p)), 
                                    rep("Asian", length(aware.p)))  
evs.c.pid7.race_south$aware.p <- aware.p  
gg.c.pid7.race_south <- ggplot(evs.c.pid7.race_south, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and Region (South)") +
  geom_text(data=subset(evs.c.pid7.race_south, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (Midwest)
c.ctrl_midwest = replace(c.ctrl, "region", "Midwest")
evs.c.pid7.race_midwest <- rbind(ans.ci(predict(reg.c.pid7_race.reg, c.ctrl_midwest, type="response")),
                                 ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_midwest, "race", "Black"), type="response")),
                                 ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_midwest, "race", "Hispanic"), type="response")),
                                 ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_midwest, "race", "Asian"), type="response")))
evs.c.pid7.race_midwest$id.group <- c(rep("White", length(aware.p)), 
                                      rep("Black", length(aware.p)), 
                                      rep("Hispanic", length(aware.p)), 
                                      rep("Asian", length(aware.p)))  
evs.c.pid7.race_midwest$aware.p <- aware.p  
gg.c.pid7.race_midwest <- ggplot(evs.c.pid7.race_midwest, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and Region (Midwest)") +
  geom_text(data=subset(evs.c.pid7.race_midwest, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (West)
c.ctrl_west = replace(c.ctrl, "region", "West")
evs.c.pid7.race_west <- rbind(ans.ci(predict(reg.c.pid7_race.reg, c.ctrl_west, type="response")),
                              ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_west, "race", "Black"), type="response")),
                              ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_west, "race", "Hispanic"), type="response")),
                              ans.ci(predict(reg.c.pid7_race.reg, replace(c.ctrl_west, "race", "Asian"), type="response")))
evs.c.pid7.race_west$id.group <- c(rep("White", length(aware.p)), 
                                   rep("Black", length(aware.p)), 
                                   rep("Hispanic", length(aware.p)), 
                                   rep("Asian", length(aware.p)))  
evs.c.pid7.race_west$aware.p <- aware.p  
gg.c.pid7.race_west <- ggplot(evs.c.pid7.race_west, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted Party Identity") + ggtitle("Race and Region (West)") +
  geom_text(data=subset(evs.c.pid7.race_west, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

c.pid.interaction_race.region = gg.c.pid7.race_north + gg.c.pid7.race_south + gg.c.pid7.race_midwest + gg.c.pid7.race_west + plot_layout(nrow = 2)

## CCES POLICY PREFERENCES ##

# CCES Party identity and interaction with Race and Region (Northeast)
c.ctrl_north = replace(c.ctrl, "region", "Northeast")
evs.c.atts.race_north <- rbind(ans.ci(predict(reg.c.atts_race.reg, c.ctrl_north, type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_north, "race", "Black"), type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_north, "race", "Hispanic"), type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_north, "race", "Asian"), type="response")))
evs.c.atts.race_north$id.group <- c(rep("White", length(aware.p)), 
                                    rep("Black", length(aware.p)), 
                                    rep("Hispanic", length(aware.p)), 
                                    rep("Asian", length(aware.p)))  
evs.c.atts.race_north$aware.p <- aware.p  
gg.c.atts.race_north <- ggplot(evs.c.atts.race_north, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Race and Region (Northeast)") +
  geom_text(data=subset(evs.c.atts.race_north, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (South)
evs.c.atts.race_south <- rbind(ans.ci(predict(reg.c.atts_race.reg, c.ctrl, type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl, "race", "Black"), type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl, "race", "Hispanic"), type="response")),
                               ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl, "race", "Asian"), type="response")))
evs.c.atts.race_south$id.group <- c(rep("White", length(aware.p)), 
                                    rep("Black", length(aware.p)), 
                                    rep("Hispanic", length(aware.p)), 
                                    rep("Asian", length(aware.p)))  
evs.c.atts.race_south$aware.p <- aware.p  
gg.c.atts.race_south <- ggplot(evs.c.atts.race_south, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Race and Region (South)") +
  geom_text(data=subset(evs.c.atts.race_south, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (Midwest)
c.ctrl_midwest = replace(c.ctrl, "region", "Midwest")
evs.c.atts.race_midwest <- rbind(ans.ci(predict(reg.c.atts_race.reg, c.ctrl_midwest, type="response")),
                                 ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_midwest, "race", "Black"), type="response")),
                                 ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_midwest, "race", "Hispanic"), type="response")),
                                 ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_midwest, "race", "Asian"), type="response")))
evs.c.atts.race_midwest$id.group <- c(rep("White", length(aware.p)), 
                                      rep("Black", length(aware.p)), 
                                      rep("Hispanic", length(aware.p)), 
                                      rep("Asian", length(aware.p)))  
evs.c.atts.race_midwest$aware.p <- aware.p  
gg.c.atts.race_midwest <- ggplot(evs.c.atts.race_midwest, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Race and Region (Midwest)") +
  geom_text(data=subset(evs.c.atts.race_midwest, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

# CCES Party identity and interaction with Race and Region (West)
c.ctrl_west = replace(c.ctrl, "region", "West")
evs.c.atts.race_west <- rbind(ans.ci(predict(reg.c.atts_race.reg, c.ctrl_west, type="response")),
                              ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_west, "race", "Black"), type="response")),
                              ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_west, "race", "Hispanic"), type="response")),
                              ans.ci(predict(reg.c.atts_race.reg, replace(c.ctrl_west, "race", "Asian"), type="response")))
evs.c.atts.race_west$id.group <- c(rep("White", length(aware.p)), 
                                   rep("Black", length(aware.p)), 
                                   rep("Hispanic", length(aware.p)), 
                                   rep("Asian", length(aware.p)))  
evs.c.atts.race_west$aware.p <- aware.p  
gg.c.atts.race_west <- ggplot(evs.c.atts.race_west, aes(x=aware.p, y=est, ymin=lo, ymax=hi, group=id.group, linetype=id.group)) +
  geom_ribbon(alpha=.2) + ylim(0,1) +
  scale_x_continuous(breaks=c(0,.25,.5,.75,1), labels=c(0,25,50,75,100)) +
  geom_line() + xlab("Political awareness percentile") + ylab("Predicted policy preferences") + ggtitle("Race and Region (West)") +
  geom_text(data=subset(evs.c.atts.race_west, aware.p==1), aes(label=id.group), hjust=0, nudge_x =0.02, lineheight=.8) +
  coord_cartesian(clip='off') +
  theme_idkn2()

c.pp.interaction_race.region = gg.c.atts.race_north + gg.c.atts.race_south + gg.c.atts.race_midwest + gg.c.atts.race_west + plot_layout(nrow = 2)
c.pp.interaction_race.region

## REGRESSION TABLE ##

stargazer(list(reg.c.atts_race.reg, reg.c.pid7_race.reg),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          column.labels = c("Policy (CCES)", "Party (CCES)"),
          dep.var.labels.include = F,
          type = "latex")

#######################################################
# RACE AND INCOME INTERACTION
#######################################################


# CCES data, party identity ** interaction race and income
reg.c.pid7_race.inc <- svyglm(
  pid7 ~
    race * aware.p +
    religion * aware.p +
    female * aware.p +
    lgbt * aware.p +
    union.member * aware.p +
    race * factor(income) * aware.p +
    veteran * aware.p +
    age + married + factor(income) + educ + rel.atd +
    region + year,
  w.cces
)

# CCES data, policy preferences ** interaction race and income
reg.c.atts_race.inc <- svyglm(
  atts ~
    race * aware.p +
    religion * aware.p +
    female * aware.p +
    lgbt * aware.p +
    union.member * aware.p +
    veteran * aware.p +
    race * factor(income) * aware.p +
    age + married + factor(income) + educ + rel.atd +
    region + year,
  w.cces
)

# ANES data, party identity ** interaction race and income
reg.a.pid7_race.inc <- svyglm(
  pid7 ~
    race * aware.p +
    religion * aware.p +
    female * aware.p +
    lgbt * aware.p +
    union.member * aware.p +
    veteran * aware.p +
    race * factor(income) * aware.p +
    age + married + factor(income) + educ + rel.atd +
    region + year,
  w.anes
)

# ANES data, policy preferences ** interaction race*income
reg.a.atts_race.inc <- svyglm(
  atts ~
    race * aware.p +
    religion * aware.p +
    female * aware.p +
    lgbt * aware.p +
    union.member * aware.p +
    veteran * aware.p +
    race * factor(income) * aware.p +
    age + married + factor(income) + educ + rel.atd +
    region + year,
  w.anes
)


### make one plot per race, with income level varied in each plot

c.ctrl_race = c.ctrl ## starts at white
c.ctrl_white = c.ctrl_race
c.ctrl_black = replace(c.ctrl_race, "race", "Black")
c.ctrl_hispanic = replace(c.ctrl_race, "race", "Hispanic")
c.ctrl_asian = replace(c.ctrl_race, "race", "Asian")


### PARTY ID
evs.c.pid7.white_income = rbind(
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_white , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_white , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_white , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_white , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_white , "income" , 5), type = "response"
  ))
)
evs.c.pid7.white_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.pid7.white_income$aware.p <- aware.p

gg.c.pid7.white_income <-
  ggplot(
    evs.c.pid7.white_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, colour = NA) +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Party Identity") + ggtitle("Race and income (White)") +
  geom_text(
    data = subset(evs.c.pid7.white_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') + 
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

white = gg.c.pid7.white_income +  guides(color = "none")

evs.c.pid7.black_income = rbind(
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_black , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_black , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_black , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_black , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_black , "income" , 5), type = "response"
  ))
)
evs.c.pid7.black_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.pid7.black_income$aware.p <- aware.p

gg.c.pid7.black_income <-
  ggplot(
    evs.c.pid7.black_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA)  +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Party Identity") + ggtitle("Race and income (Black)") +
  geom_text(
    data = subset(evs.c.pid7.black_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') +
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

black = gg.c.pid7.black_income  + guides(color = "none")



evs.c.pid7.hispanic_income = rbind(
  ans.ci(predict(
    reg.c.pid7_race.inc,
    replace(c.ctrl_hispanic , "income" , 1),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc,
    replace(c.ctrl_hispanic , "income" , 2),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc,
    replace(c.ctrl_hispanic , "income" , 3),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc,
    replace(c.ctrl_hispanic , "income" , 4),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc,
    replace(c.ctrl_hispanic , "income" , 5),
    type = "response"
  ))
)
evs.c.pid7.hispanic_income$id.group <-
  c(
    rep("One", length(aware.p)),
    rep("Two", length(aware.p)),
    rep("Three", length(aware.p)),
    rep("Four", length(aware.p)),
    rep("Five", length(aware.p))
  )
evs.c.pid7.hispanic_income$aware.p <- aware.p

gg.c.pid7.hispanic_income <-
  ggplot(
    evs.c.pid7.hispanic_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA)  +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Party Identity") + ggtitle("Race and income (Hispanic)") +
  geom_text(
    data = subset(evs.c.pid7.hispanic_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') +
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

hispanic = gg.c.pid7.hispanic_income + guides(color = "none")


evs.c.pid7.asian_income = rbind(
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_asian , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_asian , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_asian , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_asian , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.pid7_race.inc, replace(c.ctrl_asian , "income" , 5), type = "response"
  ))
)
evs.c.pid7.asian_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.pid7.asian_income$aware.p <- aware.p

gg.c.pid7.asian_income <-
  ggplot(
    evs.c.pid7.asian_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA) +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Party Identity") + ggtitle("Race and income (Asian)") +
  geom_text(
    data = subset(evs.c.pid7.asian_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') +
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")


asian = gg.c.pid7.asian_income + guides(color = "none")

c.id.interaction_race.income = (white +
                                  theme(plot.margin = unit(c(0, 30, 0, 0), "pt"))
                                + black + hispanic + asian)  +
  plot_layout(guides = "collect")


### POLICY PREFS
evs.c.atts.white_income = rbind(
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_white , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_white , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_white , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_white , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_white , "income" , 5), type = "response"
  ))
)
evs.c.atts.white_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.atts.white_income$aware.p <- aware.p

gg.c.atts.white_income <-
  ggplot(
    evs.c.atts.white_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA) +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") +
  ylab("Predicted Policy Preferences") + ggtitle("Race and income (White)") +
  geom_text(
    data = subset(evs.c.atts.white_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') + 
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

white_atts = gg.c.atts.white_income + guides(color = "none")

evs.c.atts.black_income = rbind(
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_black , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_black , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_black , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_black , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_black , "income" , 5), type = "response"
  ))
)
evs.c.atts.black_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.atts.black_income$aware.p <- aware.p

gg.c.atts.black_income <-
  ggplot(
    evs.c.atts.black_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA)  +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Policy Preferences") + ggtitle("Race and income (Black)") +
  geom_text(
    data = subset(evs.c.atts.black_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') + 
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

black_atts = gg.c.atts.black_income + guides(color = "none")



evs.c.atts.hispanic_income = rbind(
  ans.ci(predict(
    reg.c.atts_race.inc,
    replace(c.ctrl_hispanic , "income" , 1),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc,
    replace(c.ctrl_hispanic , "income" , 2),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc,
    replace(c.ctrl_hispanic , "income" , 3),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc,
    replace(c.ctrl_hispanic , "income" , 4),
    type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc,
    replace(c.ctrl_hispanic , "income" , 5),
    type = "response"
  ))
)
evs.c.atts.hispanic_income$id.group <-
  c(
    rep("One", length(aware.p)),
    rep("Two", length(aware.p)),
    rep("Three", length(aware.p)),
    rep("Four", length(aware.p)),
    rep("Five", length(aware.p))
  )
evs.c.atts.hispanic_income$aware.p <- aware.p

gg.c.atts.hispanic_income <-
  ggplot(
    evs.c.atts.hispanic_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA)  +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Policy Preferences") + ggtitle("Race and income (Hispanic)") +
  geom_text(
    data = subset(evs.c.atts.hispanic_income, aware.p == 1),
    aes(label = id.group),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') + 
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

hispanic_atts = gg.c.atts.hispanic_income + guides(color = "none")


evs.c.atts.asian_income = rbind(
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_asian , "income" , 1), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_asian , "income" , 2), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_asian , "income" , 3), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_asian , "income" , 4), type = "response"
  )),
  ans.ci(predict(
    reg.c.atts_race.inc, replace(c.ctrl_asian , "income" , 5), type = "response"
  ))
)
evs.c.atts.asian_income$id.group <- c(
  rep("One", length(aware.p)),
  rep("Two", length(aware.p)),
  rep("Three", length(aware.p)),
  rep("Four", length(aware.p)),
  rep("Five", length(aware.p))
)
evs.c.atts.asian_income$aware.p <- aware.p

gg.c.atts.asian_income <-
  ggplot(
    evs.c.atts.asian_income,
    aes(
      x = aware.p,
      y = est,
      ymin = lo,
      ymax = hi,
      group = id.group,
      color = id.group,
      fill = id.group
    )
  ) +
  geom_ribbon(alpha = .2, color = NA) +
  scale_x_continuous(breaks = c(0, .25, .5, .75, 1),
                     labels = c(0, 25, 50, 75, 100)) +
  geom_line() + xlab("Political awareness percentile") + 
  ylab("Predicted Policy Preferences") + ggtitle("Race and income (Asian)") +
  geom_text(
    data = subset(evs.c.atts.asian_income, aware.p == 1),
    aes(label = id.group,),
    hjust = 0,
    nudge_x = 0.02,
    lineheight = .8
  ) +
  coord_cartesian(clip = 'off') + 
  scale_fill_discrete(breaks = c('One', 'Two', 'Three', 'Four', 'Five'),
                      name = "Income Level")

asian_atts = gg.c.atts.asian_income + guides(color = "none")

c.pp.interaction_race.income = (white_atts + 
                                  theme(plot.margin = unit(c(0, 30, 0, 0), "pt")) + 
                                  black_atts + hispanic_atts + asian_atts) +
  plot_layout(guides = "collect")

#############
#REGRESSION TABLE PREDICTING INCOME INTERACTION
stargazer(list(reg.c.atts_race.inc, reg.c.pid7_race.inc),
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          column.labels = c("Policy (CCES)", "Party (CCES)"),
          dep.var.labels.include = F,
          type = "latex")


#####################################################
# INTRA GROUP VARIATION DENSITY PLOTS
####################################################
library(stats)
quantile(cces$aware.p, probs = c(0.25,0.5,0.75,0.9))
quantile(anes$aware.p, probs = c(0.25,0.5,0.75,0.9))

#CCES#
v.cces <- cces %>% 
  mutate(v.aware.p = 
           ifelse(cces$aware.p <= 0.2417333, "low", 
                  ifelse(cces$aware.p >= 0.2417333 & cces$aware.p <= 0.7727864, "med", 
                         ifelse(cces$aware.p >= 0.7727864, "high", NA))))

atts.var.plots.c <- v.cces %>% 
  ggplot(aes(x = atts, color=v.aware.p, weights=weight)) +
  geom_density() +
  facet_wrap(vars(race)) + 
  labs(x = "Policy Preferences",
       title = "CCES, Policy Preferences",
       caption = "Policy preferences: 0 = Very Conservative, 1 = Very Liberal") + 
  theme_bw() +
  scale_color_discrete(name  ="Political Awareness Level",
                       breaks=c("low", "med", "high"),
                       labels=c("Low", "Medium", "High"))


###########################################################
# VARIANCE FUNCTION REGRESSIONS
###########################################################
library(data.table)

## Create a linear regression
#POLICY PREFERENCES 
reg.c.atts <- svyglm(atts~
                       race*aware.p+
                       religion*aware.p+
                       female*aware.p+
                       lgbt*aware.p+
                       union.member*aware.p+
                       veteran*aware.p+ age+married+income+educ+rel.atd+region+year, w.cces, weights = weight)

#PARTY IDENTITY
reg.c.pid <- svyglm(pid7~
                      race*aware.p+
                      religion*aware.p+
                      female*aware.p+
                      lgbt*aware.p+
                      union.member*aware.p+
                      veteran*aware.p+ age+married+income+educ+rel.atd+region+year, w.cces, weights = weight)

#Extract the residuals and square them
residuals_atts <- reg.c.atts$residuals
res_sq_atts <- (residuals_atts)^2

residuals_pid <- reg.c.pid$residuals
res_sq_pid <- (residuals_pid)^2

#create new dataframes of the model results + a column for the squared residuals
cces.regression.atts <- reg.c.atts$model
cces.res.atts <- cces.regression.atts %>% 
  mutate(res_sq_atts = res_sq_atts)
setnames(cces.res.atts, old=c("(weights)"), new=c("weight"))

cces.regression.pid <- reg.c.pid$model
cces.res.pid <- cces.regression.pid %>% 
  mutate(res_sq_pid = res_sq_pid)
setnames(cces.res.pid, old=c("(weights)"), new=c("weight"))

##run these data through linear regressions with a gamma distribution
#Policy preferences
atts.g <- glm(res_sq_atts~
                race*aware.p+
                religion*aware.p+
                female*aware.p+
                lgbt*aware.p+
                union.member*aware.p+
                veteran*aware.p+
                age+married+income+educ+rel.atd+region+year, cces.res.atts, family = Gamma(link="log"), weights = weight)

#Party identity
pid.g <- glm(res_sq_pid~
               race*aware.p+
               religion*aware.p+
               female*aware.p+
               lgbt*aware.p+
               union.member*aware.p+
               veteran*aware.p+
               age+married+income+educ+rel.atd+region+year, cces.res.pid, family = Gamma(link="log"), weights = weight)

#Table
stargazer(atts.g, pid.g,
          single.row = T,
          stars = c(.001, .01, .05, .1),
          fontsize = "small",
          dep.var.labels.include = F,
          column.labels = c("Policy Preferences", "Party Identity"),
          type = "latex"
)